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Abstract 


In this work, we consider a conventional test of gravitational wave (GW) propagation which is based on the 
phenomenological parameterized dispersion relation to describe potential departures from General Relativity (GR) 
along the propagation of GWs. But different from tests conventionally performed previously, we vary multiple 
deformation coefficients simultaneously and employ the principal component analysis (PCA) method to remedy 
the strong degeneracy among deformation coefficients and obtain informative posteriors. The dominant PCA 
components can be better measured and constrained, and thus are expected to be more sensitive to potential 
departures from the waveform model. Using this method we analyze ten selected events and get the result that the 
combined posteriors of the dominant PCA parameters are consistent with GR within 99.7% credible intervals. The 
standard deviation of the first dominant PCA parameter is three times smaller than that of the original dispersion 
parameter of the leading order. However, the multi-parameter test with PCA is more sensitive to not only potential 
deviations from GR but also systematic errors of waveform models. The difference in results obtained by using 
different waveform templates hints that the demands of waveform accuracy are higher to perform the multi- 
parameter test with PCA. Whereas, it cannot be strictly proven that the deviation is indeed and only induced by 
systematic errors. It requires more thorough research in the future to exclude other possible reasons in parameter 
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1. Introduction 


General Relativity (GR) is considered as the most successful 
gravity theory which has been intensively tested in the past 
across scales of laboratory experiments to observations of the 
large-scale structure of the universe (Adelberger 2001; Hoyle 
et al. 2001; Stairs 2003; Jain & Khoury 2010; Wex 2014; 
Will 2014; Berti et al. 2015; Manchester 2015; Koyama 2016; 
Kramer 2017). However, on the theoretical side, difficulties in 
problems of singularity and quantization (DeWitt 1967; 
Kiefer 2007) hint at the possible incompleteness of GR. On 
the observational side, to explain current observations of 
galaxies and cosmology within the framework of GR, 
conceptions like dark matter and dark energy have to be 
introduced (Frieman et al. 2008; Porter et al. 2011), while there 
is another possibility that GR might be invalid at this scale 
(Debono & Smoot 2016). These facts continuously motivate 
people to pursue higher precision or develop new methods for 
performing tests on GR. In recent years, gravitational wave 
(GW) observations provide an unprecedented way to test GR. 

The direct detection of the binary black hole coalescence 
event, GW150914 (Abbott et al. 2016a), initiated the era of 
GW astronomy. More than 90 GW events from compact binary 


coalescence have been detected (Abbott et al. 2019a, 2021; 
Collaboration et al. 2021a, 2023) in the three previous 
observing runs of the LIGO, Virgo, and KAGRA collaboration 
(LVK). Thorough tests on GR have been performed by LVK 
based on these detected GW data (Collaboration et al. 
2021, 2023, 2021c). 

The tests performed by LVK include three aspects of GWs: 
generation, propagation, and polarization. The methods used 
can be classified as consistency tests and parameterized tests 
(Collaboration et al. 2021, 2023, 2021c). Consistency tests aim 
at checking whether the observed data are consistent with 
predictions of GR, such as the residual test where the best-fit 
waveform will be subtracted from data and checking whether 
there is remnant coherent power in residuals (Cornish & 
Littenberg 2015), or the inspiral-merger-ringdown consistency 
test where the low frequency and high frequency parts of 
signals are used to perform parameter estimation separately and 
checking whether the results are consistent (Ghosh et al. 
2016, 2017). In addition, parameterized tests adopt specific 
parameterization of possible deviations from GR. For example, 
the parameterized test of GW generation (Li et al. 2012; 
Agathos et al. 2014) utilizes the parameterization based on the 
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post-Newtonian (PN) expansion structure. Deformation coeffi- 
cients are added to original PN coefficients which are solely 
determined by masses and spins of the binary in GR. The 
additional deformation coefficients will be estimated as free 
parameters in Bayesian parameter estimation, which can 
capture various possible unmodeled effects and can also be 
mapped to modifications in specific alternative gravity theories 
through the parameterized post-Einsteinian (ppE) framework 
(Yunes & Pretorius 2009). Similarly, the parameterized 
dispersion relation is used in the test of GW propagation 
(Mirshekar et al. 2012). Additional power terms of GW 
momentum are added in the dispersion relation, which can 
make different frequency components of GWs propagate with 
different speeds, thus distorting the observed waveform 
(Will 1998; Mirshekari et al. 2012). 

There are multiple deformation coefficients in parameterized 
tests for capturing various potential deviations. Due to the 
limitation of sensitivity of current detectors, if varying all 
deformation coefficients simultaneously in parameter estima- 
tion, correlations among parameters will lead to uninformative 
posteriors. Thus, in tests performed by LVK, only one 
deformation parameter is allowed to vary at a time (Abbott 
et al. 2016b, 2019b; Collaboration et al. 2021, 2023, 2021c). 

However, in the most general and agnostic-priori case, all 
deformation parameters need to be considered as free 
parameters in Bayesian inference. The class of parameterized 
tests where multiple deformation parameters are constrained 
simultaneously is referred to as multi-parameter tests (Datta 
et al. 2021; Gupta et al. 2020; Saleem et al. 2022). Previous 
works (Datta et al. 2021; Gupta et al. 2020) show that the 
multi-band observations by third-generation ground-based 
detectors and space-borne detectors can allow people to 
perform multi-parameter tests and get tight constraints on PN 
deformation coefficients. Another approach to perform multi- 
parameter tests is employing the method of Principal Comp- 
onent Analysis (PCA) to reduce correlations among parameters 
and get informative posteriors (Pai & Arun 2012; Ohme et al. 
2013; Saleem et al. 2022; Datta et al. 2024; Datta 2023; Shoom 
et al. 2023). The previous work (Saleem et al. 2022) applies 
this method in the parameterized test of GW generation. In this 
work, we will extend this method to the parameterized test of 
GW propagation. The rest of this paper is organized as follows. 
In Section 2, we introduce the parameterized test of GW 
propagation, and the method of PCA for multi-parameter tests. 
We apply this method to 10 selected GW events, and present 
results in Section 3. The summary is shown in Section 4. 


2. Methods 
2.1. The Multi-parameter Test of GW Propagation 


The main effort of this paper is extending the multi- 
parameter test of GW generation with PCA performed in the 
previous work (Saleem et al. 2022) to the test of propagation. 
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We focus on the effects of dispersion, while other propagation 
effects like birefringence (Zhao et al. 2020; Wang et al. 2021b; 
Okounkova et al. 2021; Niu et al. 2022; Wang et al. 2022) or 
amplitude damping (Belgacem et al. 2018; Nishizawa 2018) 
are not included in this paper. The test of propagation 
considered here is based on the phenomenological modified 
dispersion relation (Mirshekari et al. 2012) which reads 


E? = p? + M Agp*, (1) 


a 


where E and p are the energy and momentum of GW 
respectively, A, are the free parameters for capturing various 
potential deviations from GR, and a values correspond to 
modifications at different frequency orders. When A, = 0 for 
all a, this modified dispersion relation returns to the case of 
GR. Leading-order modifications in various alternative theories 
can be mapped to the terms with different values of a (Yunes & 
Pretorius 2009; Mirshekari et al. 2012; Yunes et al. 2016). 

This phenomenological modified dispersion relation is a 
generic framework and can cover a wide variety of alternative 
theories, where the behaviors of GW propagation are different 
from GR, including theories with massive gravitons or various 
Lorentz-violating alternative gravity theories (Mirshekari et al. 
2012). For example, the case of œa = 0 and A, > 0 corresponds 
to deformations induced by massive gravitons (Will 1998); the 
leading-order corrections in multi-fractal spacetime (Cal- 
cagni 2010) and the doubly special relativity (Amelino- 
Camelia 2002) can be mapped to the term of a= 2.5 and 3; 
The Horáva-Lifshitz theory (Hofava 2009), extra-dimensional 
theories (Sefiedgar et al. 2011), and the standard model 
extension where if only non-birefringent effects are considered 
(Kostelecky & Mewes 2016), have leading modifications 
with a=4. 

Additional power terms of momentum in the modified 
dispersion relation will make different frequency components 
of GW propagate with different speeds, thus leaving imprints in 
the waveform observed. Assuming the waveform in the local 
wave zone is consistent with GR, the additional corrections in 
phase induced by dispersion effects along propagation are 
given by Mirshekari et al. (2012) 


a—1 
TUI e, os, 
a-l : 


= In (7. Mf), 


a 


69,(f) = sign(Aa) 


a= 1. 


(2) 


In the above equation, M is the chirp mass of the binary, and 
Na = h|A.|"/ (^72 with the Planck constant A, Da is defined as 


l-a z ^1a—2 
D,= Gt (Ore ow. O) 
Ho 79 QOO + zy + Op 
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where Ho, is the Hubble constant, z is the redshift of the 
source, and Qm and ©, are the matter and dark energy 
density parameters respectively. We use the values 
reported in Planck 2018 results (Aghanim et al. 2020), where 
Ho = 67.66 km s! Mpc !, Qu — 0.3111 and N4 = 0.6889. 

For the convenience of varying all deformation coefficients 
simultaneously in parameter estimation, we consider a para- 
meterization which is slightly different from that used by LVK 
(Abbott et al. 2016b, 2019b; Collaboration et al. 2021, 
2023, 2021c). Dimensionless coefficients 


p, = 1 Gpe(10? Mo)! -^A27?, (4) 


are constructed, where the factors of distance and mass in front 
are adjusted according to typical values of detected sources in 
order to scale the magnitude of deformation coefficients to a 
close order. Therefore, the phase correction of Equation (2) can 
be rewritten as 


6&,(f) = sign(Aa) 
a+! (10 Ms\*!/ n, "T 
6$, CES S (mme y (yate, o 


bb, 


r In(nMf), a — 1. 


D, 
1 Gpc 
(5) 


The above modifications are added into a GR waveform 
model as the template used in the Bayesian parameter 
estimation. The measurement of all free parameters in the 
template, including deformation coefficients of dispersion 
effects and source properties in GR, will be presented through 
posterior distributions obtained by the nested sampling 
(Skilling 2004, 2006) or the Markov chain Monte Carlo 
(MCMC) (Hastings 1970; Foreman-mackey et al. 2013) 
algorithms. The posterior samples of 6¢,, after marginalizing 
over GR parameters, are manipulated by the PCA method to 
reduce correlations among these coefficients and obtain 
informative posteriors. We will describe the above procedures 
in more detail in the following. 


2.2. Bayesian Parameter Estimation with GWs 


Next, we present a brief review of Bayesian parameter 
estimation for GW data (Abbott et al. 2020). The physical 
information is extracted from observations through the 
Bayesian framework which estimates the probability distribu- 
tions of parameters 0 in the model M according to observed 
data d. The measurement is encoded by posterior distributions 
p(0|d, M). According to Bayes' theorem, the posterior is given 
by 


p(0)p(d|0, M) 
p(d) 


Here, p(0) is the prior probability which reflects our prior 
knowledge of the parameters ahead of observations. The prior 


p(0d, M) = (6) 
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Table 1 

Priors Used in Parameter Estimation 
Parameters Priors 
ai Uniform [0, 0.99] 
a Uniform [0, 0.99] 
0 Sine [0, 7] 
05 Sine [0, 7] 
On Uniform [0, 27] 
$12 Uniform [0, 27] 
w Uniform [0, 7] 


a 
Q 


osine [—7/2, 7/2] 


a Uniform [0, 27] 

On Sine [0, 7] 

o Uniform [0, 27] 

ôpo Uniform [—10, 10] 
605.5 Uniform [—800, 800] 
663 Uniform [—800, 800] 
604 Uniform [—500, 500] 


Note. The notation follows symbols used in Bibly (Ashton et al. 2019; 
Romero-Shaw et al. 2020). For parameters of masses, luminosity distance, and 
merger time, we used uniform priors but the ranges are adjusted according to 
results published by LVK to reduce the computation burden of stochastic 
sampling. However, we believe enough large parameter spaces have been 
explored by ensuring that posteriors decline to zero at boundaries of priors. For 
waveform models of aligned spins, we use the uniform prior in the range of 
[—0.99, 0.99] for the magnitude of projection of spin onto the orbital angular 
momentum. And for the waveform model NRSur7dq4, we adjust the prior 
range of mass ratio and magnitude of spins for satisfying the requirements of 
the waveform model. 


distributions of all parameters used in the parameter estimation 
are shown in Table 1. p(d|0, M) is the likelihood whose value is 
the probability of the occurrence of a noise realization that is 
just equal to the observed data d subtracting a GW signal given 
by the model M with the specific values of 0. p(d) is the 
normalizing factor and also is referred to as evidence for its 
usage in model selection. 

In the GW data analysis context, assuming the noise is 
Gaussian and stationary, the likelihood can be written as 


uM» 


s (7) 


p(d|0, M) x exp] 
Here, h is the GW signal given by a set of parameters 0 in a 
specific model M, and the brackets denote the inner product 
which is defined as 
*(f)b 
(alb) = 4Re fap CE HEY) (8) 
S,CÉ) 
where S$,(f) is the power spectral density (PSD) of detector 
noise. For GW signals, the dimension of parameter space is 
usually high, and the likelihood evaluation is usually 
computationally expensive. Computing posteriors on a grid is 
impractical. Stochastic sampling algorithms, including nested 
sampling (Skilling 2004, 2006) or the MCMC (Hastings 1970; 
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Foreman-Mackey et al. 2013), are used to obtain random 
samples whose densities can be approximated to the posterior 
probabilities. 

In our situation, the model contains multiple dispersion 
parameters 6d, together with GR parameters. Different from 
tests performed by LVK, we vary multiple deformation 
parameters simultaneously in the parameter estimation intro- 
duced above. After marginalizing GR parameters, the posterior 
samples of dispersion parameters are manipulated by the 
procedure discussed in the next subsection to reduce 
correlations. 


2.3. Principal Component Analysis 


Different from the tests performed by LVK, multiple 
dispersion parameters are varied simultaneously in our 
parameter estimation. Due to correlations among these 
deformation parameters, varying multiple parameters will lead 
to less informative posteriors. However, previous works (Pai & 
Arun 2012; Ohme et al. 2013; Saleem et al. 2022; Datta et al. 
2024; Datta 2023; Shoom et al. 2023; Niu et al. 2024) show 
that this problem can be remedied by the method of PCA. 

The method of PCA will find a new set of bases for the 
parameter space by a linear combination of original parameters. 
The newly constructed parameters have reduced correlations, 
and can be better constrained by data. To obtain the new bases, 
we first need to compute the covariance matrix X for the 
posterior samples of dispersion parameters 6¢, after margin- 
alizing GR parameters, which is given by 


Ly = (60; — (6;)) (6d; — (66). a 


Here, (---) denotes the expectation value. Then, upon 
diagonalizing the covariance matrix X, it can be written in 
terms of eigenvectors and eigenvalues as 


= = VAU". (10) 


A is a diagonal matrix whose diagonal elements are 
eigenvalues of X. U is a matrix composed by eigenvectors of 
X. The eigenvectors are the new bases we are looking for, in 
which each dimension is linearly uncorrelated. The posteriors 
on this new set of bases can be obtained by the transformation 


6QPC^ — Y Ullóo., (11) 
j 


where sprei denotes the i-th principal component. 

The levels signifying how principal a component is are 
indicated by the values of corresponding eigenvalues. The 
eigenvector with the smallest eigenvalue is the most dominant 
component corresponding to the new PCA parameter with the 
smallest error bar. The method of PCA redistributes the 
information on posteriors among the new bases, and collects 
information on all PN orders into dominant components. 
Therefore, potential deviations can be more obvious when 
consiering posteriors with their corresponding PCA parameters. 
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2.4. Combination of Results from Multiple Events 


According to the method introduced above, we can perform 
the multi-parameter test of GW dispersion using a single GW 
event. In this subsection, we introduce the operation for 
combining information from multiple events. 

The deviations from GR, if they exist, are believed to be 
universal for all events. Assuming the observation of each 
event is an independent measurement of deformation para- 
meters óa, we can directly multiply the posterior probabilities 
of óa given by each event to get the combined posteriors, 
which is the same process as in parameterized tests performed 
by LVK (Abbott et al. 2016b, 2019b; Collaboration et al. 
2021, 2023, 2021c). However, since the linear combination of 
original dispersion parameters to get PCA parameters sprei is 
unique for each event, it is inappropriate to multiply posteriors 
of PCA parameters. Therefore, we resample the combined 
posteriors of 6¢, and perform a separate operation of PCA to 
obtain the combined posteriors of 6QPCA. 

In practice, we use Gaussian kernel density estimation 
(KDE) to obtain posterior probabilities of ôa from its posterior 
samples for each event. Then, using a new MCMC sampling to 
sample the cumulative production of KDEs of all events, we 
can obtain the posterior samples of 6¢, combining information 
from all events. The combined posteriors of PCA parameters 
ip are obtained by performing a separate operation of PCA 
on the combined posterior samples of ġa 

As will be discussed in the next section, to verify the 
robustness of using different waveform models, we do the same 
parameter estimation with different templates for each event. 
We rely on the method discussed in Ashton & Khan (2019) to 
average the results obtained by using different waveform 
approximations, where the posterior for a single event 
estimated by using a set of models M; is given by 


p(8. (M3) = Y: pX8. M)£, (2) 


and the weights £; are defined as €;= Z;/? ;Z; with Z; denoting 
the evidence of each waveform model. 


3. Results and Discussions 


We apply the method elaborated on in the last section to 
current observed GW data and present results in this section. 
Due to the limitation on sensitivity of current detectors, 
following the multi-parameter test of generation performed in 
Saleem et al. (2022), we also do not consider the full set of 
dispersion parameters. We only include the cases of a = 0, 2.5, 
3, and 4 in our analyses by considering the examples of leading 
order modifications in alternative theories mentioned in 
Mirshekari et al. (2012); Abbott et al. (2019c). 

We analyze the 10 farthest events detected by LVK in three 
previous runs with criteria that the false alarm rate (FAR) is less 
than 107° yr, the signal-to-noise ratio (S /N) greater than 12, 
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Table 2 
The Information on 10 Selected Events Considered in Our Analyses 


Name D; (Mpc) S/N FAR (yr?) Paso 
GW170823 1940900 12.2702 «107 1.00 
GW190408 181802 1540230 14.6703 «10^? 1.00 
GW190503 185404 1520+630 122704 «10? 1.00 
GW190512 180714 14607255 12,7397 «10? 1.00 
GW190513 205428 22107910 12.5*03 1.3 x 10? 1.00 
GW190519 153544 260071729 15.9793 «105 1.00 
GW190602 175927 284071533 13.2702 «10^ 1.00 
GW190706 222641 363072000 13402 5.0 x 10? 1.00 
GW190828_063405 20707839 16.5703 «10^? 1.00 
GW190915 235702 17507248 13.1*03 «10^ 1.00 


Note. These events are 10 farthest events with criteria of FAR < 107° yr |, S/ 
N> 12, and Pasto 0.99 in the Gravitational Wave Transient Catalog 
(GWTC). The values come from gwosc.org which may be different from the 
values shown in publications of LVK (Abbott et al. 2019a, 2021; Collaboration 
et al. 2021a, 2023). 


and the probability of astrophysical origin (Pastro) greater than 
0.99. The potential dephasing caused by dispersion will 
accumulate along the propagation of GW. As can be seen 
through Equation (5), if the sensitivity to waveform dephasing 
6®, is fixed, we can generally constrain 6¢, better by sources 
with farther distance. However, this is not strictly valid in any 
cases. The correlations among parameters are complicated, and 
the sensitivity to dephasing will depend on various factors, like 
S/N, and others. Only qualitatively can we expect that better 
constraints be obtained by further sources. The information on 
selected events is summarized in Table 2. We obtain strain data 
and noise PSD data from public data release of LVK.’ 
The strain data with glitch subtraction are used in 
parameter estimation for events GW190503_185404 and 
GW190513_ 205428. 

We use bilby (Ashton et al. 2019) with nested sampler 
pymultinest (Feroz & Hobson 2007; Feroz et al. 2008; 
Buchner et al. 2014) to perform Bayesian parameter estimation. 
Since the method of PCA is more sensitive to not only possible 
deviations from GR but also systematic errors of waveform 
template, we also use different waveforms, including 
SEOBNRv4_ROM (Pürrer 2014, 2015), IMRPhenomXAS 
(Pratten et al. 2020), IMRPhenomXHM (Pratten et al. 2021), 
and NRSur7dq4 (Varma et al. 2019), in our parameter 
estimation. 

The final results are presented in Figure 1 where the 
posteriors are obtained by averaging results of different 
waveform models with weights of evidence and combining 
results of all selected events using the method discussed in 
Section 2.4, and the results of each event are also shown in 


*  https://gwosc.org/eventapi/html/GWTC/ 
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Figure 4. The posteriors of original dispersion parameters óa 
are shown on the left sides with blue shadow in each subplot, 
and posteriors of PCA parameters 69 PA are shown on the right 
sides with orange shadow. The blue and orange solid lines 
signify the 99.7% credible intervals, and the red dashed lines 
mean the GR values. Due to the limitation of plot ranges, the 
error bars of G5, 6$» 5, 663, and the upper limit of ġo are 
out of the visible ranges. The explicit values of error bars and 
departures of maximum likelihood values from GR values are 
summarized in Table 3 for quantitative comparison. 

We can observe that the posteriors of all parameters, 
including original dispersion parameters and PCA parameters, 
are consistent with GR within 99.7% credible intervals through 
Figure | and Table 3. As discussed in Section 2.3, the dominant 
PCA parameters can be better constrained, thus having smaller 
error bars, which are obviously shown in Figure | by observing 
the first three components. For example, the range of 99.7% 
credible intervals of the first dominant PCA parameter is 3 
times smaller than that of the original dispersion parameter of 
the leading order. We can also notice that the maximum 
likelihood values are not exactly the GR values, which is 
reasonable considering unavoidable errors in waveform models 
and non-perfect Gaussian noise in reality. These departures are 
more obvious in dominant PCA parameters. Observing the 
subplot of the first dominant component M and the leading 
order dispersion parameters ġo, the departure between the 
maximum likelihood value and the GR value of aan is 
2.850, while this value is 1.340 for ġo. This on one hand 
shows that the dominant PCA parameter is more sensitive to 
any departures in posteriors, and on the other hand indicates the 
demand of waveform accuracy is higher when searching for 
deviations using PCA. 

The method of PCA is more sensitive to any departures from 
the considered waveform model, including both potential 
deviations from GR and systematic errors of the template or 
unmodeled effects. Therefore, we use different waveform 
models in parameter estimation and the obtained results are 
shown in Figures 2 and 3 for posteriors of original dispersion 
parameters and PCA parameters. From these two corner plots, 
we can find that, as introduced in Section 2.3, the new bases 
constructed by PCA have less correlations. Additionally, we 
can also observe that the results given by IMRPhenomXPHM, 
which is a phenomenological waveform model incorporating 
effects of higher modes and precession induced by in-plane 
spins, have slight differences from results produced by other 
waveform templates in terms of posteriors of original 
dispersion parameters. But after the operation of PCA, the 
result of IMRPhenomXPHM has obvious differences with 
others especially in the first dominant component. This in 
another aspect shows that the operation of PCA is sensitive to 
minor changes in posteriors of original parameters, which 
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Figure 1. Posteriors of deformation parameters before and after performing PCA. The result shown here is obtained by combining posteriors of all selected events and 
averaging with weights of evidence for different waveform templates. The posteriors of original dispersion parameters are shown on the left side with blue shadows, 
and posteriors of PCA parameters are shown on the right side with orange shadows. The solid lines signify 99.7% credible intervals and the red dashed lines mean GR 


values. 
Table 3 
99.7% Credible Intervals and Distance between GR Values and Maximum Likelihood Values of Combined Results for Original Dispersion Parameters and PCA 
Parameters 

Original Dispersion Parameters PCA Parameters 

99.7% Credible Interval GR Value At 99.7% Credible Interval GR Value At 
boo [—0.154, 0.363] -1.340 Oh [—0.159 0.005 20] 2.850 
6» s [—65.1, 57.9] —0.01120 borca [—0.334 0.466] —0.480c 
65 [—41.2, 50.9] —0.06620 505A, [—2.80 5.06] —0.889c 
694 [-9.15, 6.77] 0.1900 695€, [771.3 83.0] —0.0181c 


Note. The dominant PCA parameters can be better measured and constrained, thus having smaller error bars. The operation of PCA will make the behavior of 
deviating from zero in posteriors more obvious by transforming posterior samples into the set of new bases. 


increases the demand for waveform accuracy when implement- 
ing the method of PCA. 

However, the most common reason that induces the 
deviation is the insufficient accuracy of the waveform model, 
and the different results given by different waveform models 
also hint that. Other researches also encounter similar 
deviations and attribute the deviations to possible systematic 
errors. For example, see the discussion in Appendix C of 
Collaboration et al. (2021), and Appendix A of Perkins et al. 


(2021). It cannot be strictly proven that the deviation is indeed 
and only induced by systematic errors of waveform models. 
For example, the inclusion of effects of higher modes and 
precession in the waveform model IMRPhenomXPHM may 
help to break degeneracy among parameters which can allow 
us to get better measurements of source parameters. While this 
also increases the sensitivity of non-stationary and non- 
Gaussian features in the noise realization which are also 
potential reasons to induce bias in parameter estimation, it 
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Figure 2. Posteriors of original dispersion parameters, which are combined results of all selected events, obtained by using different waveform models. The orange 


solid lines signify GR values and the dashed lines mean 99.7% credible intervals. 


requires more thorough research in the future to exclude other 
possible reasons in parameter estimation and data processing. 


4. Summary 


GW observations have opened a new era of exploring the 
nature of the universe. Detected GW events have been leading 
paradigm-shifting in research of astrophysics, cosmology, and 
gravity. As the launch of the fourth observing run of LVK, it is 
expected that more GW events will be detected and more 
amazing discoveries may be revealed. GW observations 
provide a new approach to test gravity theories. Testing gravity 


with GWs has begun to blossom since the first GW event 
detection in 2015. Using detected events, LVK has performed 
exhaustive tests of GR (Collaboration et al. 2021, 2023, 
2021c). Except for the theory-agnostic tests performed by 
LVK, theory-specific tests also have been performed by many 
independent research groups, such as (Jana & Mohanty 2019; 
Ramos & Barausse 2019; Wang et al. 2021a; Perkins et al. 
2021; Gong et al. 2022; Wu et al. 2022; Haegel et al. 2023). 

In this work, we consider the test of the GW dispersion 
relation which is conventionally performed by LVK. But in 
tests of LVK, to avoid the correlations among deformation 
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Figure 3. The results of applying PCA on the posteriors shown in Figure 2. Same as Figure 2, the orange solid lines signify GR values and dashed lines mean 99.7% 
credible intervals. Comparing with Figure 2, we can find that correlations among these parameters are reduced, and the first three dominant components are better 
constrained. We can also observe that the slight difference in posteriors of original dispersion parameters can lead to significant differences in results after applying 
PCA, especially in the first component, which hints that the multi-parameter test with PCA may require higher accuracy of waveform templates. Although the most 
common reason to induce the deviation is the insufficient accuracy of the waveform model, as discussed in Appendix C of Collaboration et al. (2021) and Appendix A 
of Perkins et al. (2021). It cannot be strictly proven that the deviation is indeed and only induced by systematic errors of waveform models. This requires more 
thorough research in the future to exclude other possible reasons in parameter estimation and data processing. 


parameters to yield uninformative posteriors, only one 
deformation parameter is varied at a time. This is equivalent 
to considering the prior of the 6 function for the fixed 
parameters. If considering the most agnostic prior, all 
deformation parameters need to be estimated simultaneously. 
Previous works (Pai & Arun 2012; Ohme et al. 2013; Saleem 
et al. 2022; Datta et al. 2024; Datta 2023; Shoom et al. 2023) 
show that the correlations can be reduced by the method 
of PCA. Using this method, the authors considered the 


parameterized test of PN structure where multiple deformation 
parameters are allowed to vary simultaneously, but still obtain 
informative posteriors meanwhile. We extend this method to 
the test of GW dispersion relation in this work. 

The biggest difference in our analyses comparing with the 
tests of LVK is that the multiple dispersion parameters are 
varied simultaneously in Bayesian parameter estimation. Then 
we transform the obtained posterior samples into a set of new 
bases constructed by the PCA. The new PCA parameters have 
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Figure 4. Posterior distributions of original dispersion parameters for each selected event. GR values are signified by red dashed lines, and 99.746 credible intervals are 
indicated by blue solid lines. The results shown here have been averaged with the weights of evidence of different waveform models. 
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less correlations, and dominant PCA parameters can be better 
measured and constrained, and thus are expected to be more 
sensitive to potential deviations from GR. We consider 10 
farthest events selected by criteria of FAR < 10 ? yr, 
S/N > 12, and Pastro > 0.99. We consider four different GR 
waveforms and add the modifications of the dispersion effect 
into them as the template in Bayesian parameter estimation. 
The results we obtain are consistent with GR within 99.7% 
credible intervals in all cases. However, we observe that the 
maximum likelihood value of the first PCA parameter has a 
relatively large departure from the GR value. This on the one 
hand shows the dominant PCA parameter is more sensitive to 
deviations from zero in posteriors, and on the other hand hints 
that the demand of waveform accuracy is higher when using 
PCA. We also find that minor differences in posteriors of 
original parameters can lead to obvious differences in results 
after the operation of PCA, which from another aspect indicates 
multi-parameter tests with PCA require more accurate wave- 
form models. Although the systematic errors of waveform 
models are the most common reason to induce deviations, it 
cannot be strictly proven that the deviation is indeed and only 
induced by systematic errors. The non-stationary and non- 
Gaussian features in the noise realization are also potential 
reasons to induce bias in parameter estimation. This requires 
more thorough research in the future to exclude other possible 
reasons in parameter estimation and data processing. 
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